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ABSTRACT 

Under  a  non-degeneracy  assumption  on  the  non-iero  entries  of  a  kiven  sparse 
matrix,  a  polynomiaUy-bonnded  algorithm  is  presented  that  perform!  row  opera¬ 
tions  on  the  given  matrix  which  reduce  it  to  a  sparsest  possible  matrix  with  the 
same  row  space.  For  each  row  of  the  matrix,  the  algorithm  performs  samaximum 
cardinality  matching  on  the  bipartite  graph  associated  with  a  submawix  which 
is  induced  by  that  row.  The  dual  of  the  optimal  matching  then  specifies  the 
row  operations  that  will  be  performed  on  that  row.  We  else  describe  a  variant 
algorithm  that  processes  the  matrix  in  place,  thus  conserving  storage  and  time. 
The  modifications  needed  to  apply  the  algorithm  to  matrices  that  do  not  neces¬ 
sarily  satisfy  the  non-degeneracy  assumptior  dso  described.  A  particularly 
promising  application  of  this  algorithm  is  in  Auction  of  linear  constraint 
matrices. 
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1.  Introduction 

An  important  factor  in  our  present  ability  to  solve  many  large-scale  numeri¬ 
cal  problems  is  the  recognition  that  these  problems  are  nearly  always  sparse,  and 
that  taking  advantage  of  sparsity  can  turn  a  hitherto  practically  unsolvable  prob¬ 
lem  into  a  solvable  one.  Perhaps  the  best  example  of  this  is  in  large-scale  linear 
programming,  where  highly  refined  sparse  matrix  factorisation  routines  have  al¬ 
lowed  problems  with  huge  coefficient  matrices  to  be  solved  (see  e.f .,  Duff  (1080) 
or  Bunch  and  Rose  (1976)).  However,  although  sparsity  is  known  to  be  helpful, 
relatively  little  attention  seems  to  have  been  paid  to  techniques  that  economically 
increase  sparsity  (decrease  density),  thereby  improving  the  efficiency  of  sparse 
algorithms.  In  this  context,  this  paper  considers  the  Sparseness  Problem  (SP): 

Given  a  large,  sparse  system  of  linear  equations 

A *  =  6,  (1) 

find  an  equivalent  system 

Az  =  f  (2) 

which  has  the  minimum  possible  number  of  non-iero  entries  in  X 

Constraints  of  the  form  (1)  are  among  the  most  common  in  large  scale 
optimisation,  so  that  it  is  potentially  very  useful  to  solve  (SP).  Under  a  non- 
degeneracy  assumption,  we  shall  present  an  efficient  algorithm  that  solves  (SP) 
using  maximum  cardinality  matching.  Sections  2-4  will  assume  familiarity  with 
notions  of  graph  theory  and  maximum  cardinality  bipartite  matching  (see,  e. 
Lawler  (1076)).  Section  2  develops  most  of  the  machinery  needed  for  the  proof, 
and  uses  it  to  derive  an  algorithm  that  solves  a  subproblem  of  (SP).  In  Section  8 
we  use  the  algorithm  of  Section  2  to  construct  the  full  algorithm,  and  prove  that 
it  solves  (SP).  We  then  give  a  variant  algorithm  that  uses  less  space  and  show 
that  it  also  solves  (SP).  Section  4  discusses  the  modifications  necessary  to  apply 
the  algorithm  on  matrices  that  do  not  necessarily  satisfy  the  non-degeneracy 
assumption.  Finally,  Section  5  considers  further  questions  raised  by  this  research. 

2.  Transfer  ms  and  the  One  Raw  AlgeHtlnu 

In  this  paper  we  shall  assume  that  the  matrix  A  in  (1)  has  full  row  rank. 
We  know  from  linear  algebra  that  (2)  is  equivalent  to  (1)  if  and  only  if  X  =  TA 
and  ?  ss  71  for  some  square  non-singular  matrix  7.  Wh  are  aiming  for  a  general 
algorithm  that  makes  no  assumptions  about  any  special  structure  in  A,  and 
thus  can  find  7  almost  solely  from  the  spanitf  potttm  of  A  (the  positions  of 
the  non-seros  in  A).  What  can  go  wrong  in  this  aim  is  that  we  can  encounter 
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"unexpected*  cancellation.  To  illustrate,  consider  the  following  two  A’s,  with 
the  same  sparsity  pattern,  treated  with  the  same  T: 

(X  -1  1\/1  1  0  0  0\  fl  0  0  0  0\ 

**-'■1°  i  0  (o  1  1  i  i)  =  lo  1  1  1  ll; 

Vo  o  1/V0  0111/  Vo  o  l  l  1/ 

fi  -l  iyi  i  o  o  o\  fi  o  o  -l  2\ 

TA&  —  1 0  1  OllO  1  1  2  3 )  =  1 0  1  1  2  31. 

Vo  0  1/Vo  0  111/  VO  0  1  1  1/ 

In  both  cases  T  represents  the  unique  linear  transformation  that  adds  the  mul¬ 
tiples  of  rows  2  and  3  to  row  1  which  makes  sero  and  avoids  fill-in  in  <*13. 
In  the  first  case  the  sparsity  increased,  in  the  second  case  it  decreased.  The 
difficulty  is  that  A\  has  some  dependent  submatrices  that  are  not  apparent  from 
the  sparsity  pattern  alone.  The  possibility  of  this  phenomonon  makes  solving 
(SP)  too  difficult  in  general,  as  shown  by  the  following  result. 

Theorem  1*  (Stockmeyer  (1082))  (SP)  is  NP-Hard  in  general.  (See  the  Appendix 
for  the  proof.)  □ 

Thus,  to  get  a  polynomial  algorithm  for  (SP),  we  must  make  some  assump¬ 
tion  about  A.  Suppose  that  A  is  si  X  »,  and  let  R  £  {l,...,m},  C  C 
{ 1, . . .,»}.  We  denote  the  submatrix  of  A  indexed  by  the  rows  in  R  and  the 
columns  in  C  by  Arc-  The  sparsity  pattern  of  Arc  naturally  induces  a  bipartite 
graph  §rc  =  (R,C,E)  where  E  =  {(<,/)  €  R  X  C  \  a#  /  0>.  Let  M{Q)  be 
the  number  of  edges  in  a  maximum  cardinality  matching  in  the  bipartite  graph 
§.  V  |A|  as  \C\,  then  the  usual  expansion  of  det  Arc  has  at  least  one  non-sero 
term  precisely  when  M($rc)  —  |R|,  and  when  A  is  “general",  we  expect  the 
converse  of  this  to  be  true  as  well.  This  reasoning  leads  us  to  assume  henceforth 
that  A  has  the 

Matching  Property  (MP):  rank  Arc  =  M($rc)  for  all  if  and  C. 

For  example,  At  above  does  not  satisfy  (MP)  whereas  M  does.  Note  that  if  the 
entries  of  A  are  independent  algebraic  indeterminates  then  (MP)  is  satisfied. 

Since  T  must  be  non-singular,  ${T)  has  a  perfect  matching  which  we  can 
assume  without  loss  of  generality  is  the  main  diagonal.  We  can  further  assume 
that  tsi  =  1,  i  =  1,2,. ..,m  by  scaling  the  rows  of  T,  so  that  the  non-sero 
entries  in  row  i  of  T  indicate  the  multipliers  for  the  rows  to  be  added  to  row  i 
of  A.  Viewed  in  this  way,  (SP)  breaks  down  into  m  one  row  sparsity  problems 
(OR8P<),  i  a*  1, 2, . . . ,  m.  (ORSPf)  is  the  problem: 

Find  {X»,  *  yi  <}  so  that 

—  A*,*  +  X*A*,# 


(3) 
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has  the  minimum  possible  number  of  non-ieros. 

Not  all  solutions  to  (3)  are  equally  good.  Since  -we  expect  that  the  amount  of 
arithemetic  needed  to  do  the  calculations  in  (3)  depends  upon  the  number  of 
rows  with  non-sero  multipliers,  ideally  we  would  also  like  to  solve  the  8trong 
(ORSPtf): 

Among  all  optimal  solutions  to  (3),  find  one  that  minimises 

K*|x»p*o}|. 

It  is  not  clear  at  this  point  that  we  can  solve  (SP)  by  successively  solving 
(ORSP,)  for  i  =  1, 2, . . m;  nevertheless  we  shall  concentrate  on  (ORSPj)  in 
the  remainder  of  this  Section. 

A  set  of  multipliers  {  X*  |  *  >  1 }  for  (3)  when  t  =  1  defines  the  following 
index  subsets: 

tf  =  {*>l|X*^0}, 

H  =  { j  |  Jt\j  =  0  and  axj  ^  0>, 

S  =  {y  j  Jt\j  —  0  and  oiy  =  0  and  a*y  ^  0  for  some  k  €  U }, 
G=*H\JS, 

F  —  U  I  *iy  ^  0  and  au  =  0>, 

P  —  FUS  =  {j\aij~0  and  a*/  ^  0  for  some  i€lT>,  and 
Z  =  {j  |aiy  =  0>. 

That  is,  U  is  the  set  of  used  rows;  H,  the  set  of  Ut  columns,  where  a  non- zero 
was  changed  to  a  zero;  S,  the  set  of  saved  columns,  where  a  sero  that  we  would 
have  expected  to  be  filled-in  (since  a*y  ^  0)  was  not  filled-in;  G,  the  set  of  good 
columns,  where  the  entry  was  actively  manipulated  for  the  better;  F  is  the  set 
of  fliled-ln  columns;  P  is  the  set  of  potential  fill-in  columns;  and  Z  is  the  set  of 
sero  columns.  The  net  decrease  in  non-seros  in  row  1  is  then  |H|  —  |F|,  which 
we  want  to  maximise  to  solve  (ORSPi).  The  next  theorem  states  the  intuitive 
result  that  if  k  columns  are  affected  for  the  good,  then  at  least  k  independent 
rows  must  have  been  used  (we  omit  the  technical  proof). 

Theorem  3.  For  any  set  of  multipliers,  M(Qua)  —  |G|»  and  hence  rank  Aero  = 

□ 

Theorem  2  implies  in  particular  that  |E7|  >  |G|  always  holds.  17  |I7|  > 
|G|,  we  can  select  a  |G|-subset  of  U  which  perfectly  matches  to  G  and  use  the 
corresponding  square  non-singular  (by  (MP))  submatrix  of  A  to  sero  out  Aio* 
thus  achieving  the  same  result  with  lets  work.  Conversely,  if  Arc  is  a  square 
submatrix  with  a  perfect  matching,  Theorem  2  ensures  that  if  we  use  Arc  to 
sero  out  Aic,  then  G  =  C,  t.e.,  only  non-ieros  in  C  are  hit,  and  fill-in  occurs  in 
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every  position  where  it  would  be  expected.  That  is,  Theorem  2  shows  the  crucial 
fact  that  (MP)  implies  that  there  is  no  "unexpected  cancellation." 

Hence,  we  can  assume  that  the  canonical  situation  is  that  \U\  =  |G|  and 
Qua  has  a  perfect  matching.  Then  Aug  is  non-singular  by  (MP),  so  the  {  X*  } 
will  be  uniquely  determined  by  solving 

\tAug  =  M  <?•  (4) 

Equation  (4)  allows  us  to  think  of  the  {  X*  }  as  coming  from  U  and  G  rather  than 
vice  versa,  thereby  reducing  (SP)  to  the  more  combinatorial  problem  of  finding 
optimal  U  and  G. 

Thus  we  need  only  consider  all  possible  U,  and  for  each  U  consider  only  the 
G  which  match  perfectly  into  U.  There  are  potentially  many  possible  ways  to 
select  G  C  { 1, 2, . . . ,  n}  so  that  G  perfectly  matches  to  U.  The  next  theorem 
shows  that  for  a  given  U  it  suffices  to  check  only  one  such  G. 

Theorem  3.  Let  Gi  and  G2  be  two  sets  of  columns  that  perfectly  match  into 
U,  and  denote  the  set  of  hit  columns  corresponding  to  Gt-  by  Hi,  i  =  1, 2,  etc. 
Then 

Proof.  Note  that  P  depends  only  on  U  and  not  on  G{.  Then  it  is  easy  to  see 
that  =  \U\  -  |Si|  and  |F,|  =  \P\  -  |4|,  so  that  \H{\  -  |F,-|  =  |tf|  -  |P|, 
<  =  1,2.  □ 

If  we  fix  a  full-rank  matching  M  in  then  any  row  subset  U  induces  a 
unique  matched  column  subset  G  relative  to  M.  Any  such  (17,  G)  pair  will  have 
a  perfect  matching,  namely  M  restricted  to  Aug,  so  (MP)  ensures  that  Aug  will 
be  non-singular.  Hence  the  multipliers  can  be  found  as  in  (4).  Theorem  3  ensures 
that  the  best  (U,  G)  pair  from  among  this  restricted  class  of  such  pairs  will  solve 
(ORSPO. 

Letting  the  dependence  of  P  on  U  be  explicit,  through  (MP),  Theorem  2  and 
Theorem  3  we  have  reduced  the  apparently  algebraic  problem  (ORSPi)  into  the 
purely  combinatorial  one  of  maximizing  \U\  —  |P(E/)|  over  all  If  C  {  2, . . . ,  m }. 
Define  R  =  {  2, . . . ,  m  }  and  U  =  R\U.  Then 

max(|Cf|  -  |P(t/)|)  =  (m  —  1)  —  min(|P(Ef)|  +  |ZT|).  (5) 

By  definition  of  T7  and  P(t7),  every  non-sero  in  Arz  (the  sero-seetion  of  row  1 
of  A)  is  contained  in  either  a  row  in  IT  or  a  column  in  P{U).  If  we  call  rows  and 
columns  lines,  then  in  this  situation  we  say  that  Arz  is  covered  by  the  lines  in 
JJ\jP{U).  Clearly  any  such  covering  of  Arz  by  lines  can  be  written  as  V\jP{U) 
for  some  U  C  R,  so  by  (5),  finding  maxcr(|l7|  —  |P(Ef)l)  is  equivalent  to  finding 
a  minimum  covering  of  Arz  by  lines.  But  by  the  classic  theorem  of  Konig  and 
Egenrary  (see  Lawler  (1976)  p.  190),  such  a  minimum  cover  can  be  computed 
through  a  maximum  matching  in  Qrz' 
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Theorem  4.  M(Qrz)  =  miny(|P(£7)|  +  |I7|),  and  a  maximum  matching  and  a 
minimum  covering  by  lines  are  dual  combinatorial  objects.  □ 

By  the  duality  theory  of  matching  algorithms,  if  we  find  a  maximum  match¬ 
ing  in  Qrz  through  a  labelling  algorithm,  then  an  optimum  U  for  (ORSPi)  is 
the  set  of  rows  reachable  via  an  alternating  path  from  an  unmatched  row.  That 
is,  Theorem  4  shows  that  this  U  solves  the  right-hand  side  of  (5),  so  it  also  solves 
the  left-hand  side,  and  so  must  solve  (ORSPi). 

Even  better,  the  next  theorem  shows  that  the  optimum  U  defined  above 
also  solves  the  Strong  (ORSPi).  For  a  network  flow  problem  with  source  s 
and  optimal  flow  /,  define  the  standard  minimum  cut  K*  =  {»  |  there 
is  an  augmenting  path  a  — ►  t  under  /}.  Note  that  the  optimal  U  defined 
above  is  a  standard  minimum  cut  for  the  usual  way  of  solving  a  maximum 
cardinality  bipartite  matching  problem  by  converting  it  to  an  equivalent  network 
flow  problem. 

Theorem  5.  In  a  given  network,  the  standard  min  cut  is  a  subset  of  every  min 
cut.  Thus  the  standard  min  cut  is  the  same  for  every  optimal  flow,  hence  it 
is  well-defined  and  has  minimum  cardinality  among  all  min  cuts  (see  Ford  and 
Fulkerson  (1962)  p.  13  for  a  proof).  □ 

Theorems  4  and  5  together  imply  that  we  can  solve  the  Strong  (ORSPi) 
through  maximum  matching,  and  that  the  optimal  U  is  unique. 

8.  Two  Algorithms  for  (SP) 

Once  we  have  found  the  optimal  U  for  each  row  i  (say,  Ut)  through  matching, 
as  noted  above  we  can  easily  generate  the  sets  of  columns  by  choosing  G{  to 
be  the  set  of  columns  that  matches  into  If*  under  the  fixed  matching  At.  These 
(U*,  Gi)  pairs  completely  determine  the  non-sero  off-diagonal  entries  of  a  matrix 
T*  as  defined  by  (4).  The  question  arises:  is  T*  non-singular? 

To  answer  this  question,  it  is  necessary  to  investigate  what  the  uniqueness 
properties  of  the  £/,•  imply  for  the  structure  of  T* .  Define  a  directed  graph  D 
with  vertices  V  =  { 1, . . . ,  m  }  and  edges  E  =  { (Jfc,  s)  |  k  €  U  J  };  thus  D  repre¬ 
sents  the  sparsity  pattern  of  T* .  If  the  row  indices  of  A  and  T*  are  ordered 
consistent  with  the  strong  component  decomposition  of  V,  then  the  decomposi¬ 
tion  induces  a  block  lower-triangular  structure  on  T* ,  where  the  diagonal  blocks 
of  T*  correspond  to  the  strong  components  of  V. 

Theorem  6.  If  l  €  U*k  and  Jk  6  U J  then  l  €  l/J. 

Proof.  For  ease  of  notation,  let  If  =  EfJ,  XT  =  \  {<}  \  U\,  P  = 

P(Ui),  and  7  =  Z{  \  P{u\).  Thus  U  and  XT  partition  the  rows  of  the  sero- 
section  of  row  i  of  A,  and  P  and  P  partition  the  columns.  Recall  that  the  rows 
in  XT  and  the  columns  in  P  are  a  minimum  cover  of  the  sero-section  of  row  t*  of 
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A  by  lines.  Thus  A#p  =  0,  and,  since  k  €  U,  Ah7  =  0.  By  the  minimality 
of  this  cover,  the  submatrix  Ajjp  has  a  row- perfect  matching,  and  so  |ZT|  lines 
are  necessary  to  coyer  it.  Let  L*  be  the  standard  minimum  set  of  lines  covering 
the  jk*k  sero- section.  8ince  Ak-p  =  0  and  k  %U,  the  submatrix  Ajj^  is  part  of 
the  sero- section  of  row  k  and  so  must  be  covered  by  the  lines  in  L*.  Consider 
the  set  of  lines  L  =  L*  [jV  \  P.  Since  the  only  non-seros  in  the  columns  in  F 
of  the  sero- section  occur  in  rows  in  TT,  L  is  a  cover  for  the  sero-section  for 

row  k.  The  only  change  in  lines  between  L*  and  L  is  in  lines  passing  through 
Ajpp\  since  L  has  only  |ZT|  lines  passing  through  Ayp,  the  minimum  possible 
number,  L  must  also  be  a  minimum  cover.  Finally,  L  contains  at  least  as  many 
rows  as  L* ,  so  that  the  U  associated  with  L  has  at  most  as  many  rows  as  the  U 
associated  with  L*,  namely  U*k.  But  U*k  has  the  minimum  possible  number  of 
rows  for  any  minimum  cover  of  the  sero-section  of  row  k,  so  L  =»  L*.  But  this 
implies  that  Uk  C  U*.  □ 

The  conclusion  of  Theorem  6  is  precisely  that  the  graph  V  is  transitively 
closed.  This  implies  that  the  blocks  of  the  block  lower-triangular  partition  of 
T*  are  either  completely  dense  or  all  sero.  In  particular,  17  J  U  { » }  =  U*k  U  {  k  } 
for  t  and  k  in  the  same  strong  component  of  P.  These  observations  allow  us  to 
prove  the  following. 

Theorem  7.  T*  is  non-singular. 

Proof.  Since  T*  is  block  lower-triangular,  it  suffices  to  show  that  the  diagonal 
blocks  of  7*  are  non-singular.  A  typical  diagonal  block  is  indexed  by  the  vertices 
in  some  strong  component,  say  B.  As  shown  above,  the  set  B*  =  is  the 

same  for  all  •  €  B,  and  B  Q  B*.  Assume  for  convenience  that  the  fixed  matching 
X  is  such  that  row  i  matches  to  column  s',  i  =  1, . . . ,  m.  V  B  =  B*,  then  the 
diagonal  block  associated  with  B  is  clearly  just  a  re-scaling  of  (Abb)-1  (Abb 
is  non-singular  by  (MP)),  and  so  is  non-singular.  Otherwise,  let  L  —  B*  \  B. 
Then  the  diagonal  Mock  associated  with  £  is  a  re- scaling  of  the  bottom  right 
corner  of  the  matrix 

(All  _  ( All 

Abl  Abb )  \&bl  Zss/ 


(this  matrix  is  non-singular  by  (MP)).  But  it  is  well-known  that  Abb  is  non¬ 
singular  if  and  only  if  All  is  non-singular.  But  All  is  indeed  non- singular  by 
(MP).  □ 

Since  T*  is  non-singular,  we  can  use  it  to  transform  A  into  X  This  way  of 
generating  A  processes  each  row  in  parallel,  t'.e.  each  row  is  solved  relative  to  the 
original  matrix  rather  than  relative  to  a  partially  transformed  matrix.  We  call 
this  procedure  the  Parallel  Algorithm  (PA). 


Proof.  Each  row  of  A  is  made  as  sparse  as 


consider  the  Sequential 


Algorithm  (SA): 

Given  A. 

For  4  =  1,2,  ...,m  do 


j  5| 

W  yv  •  f  » 

-  1  A  -/  4  6  *  T 

|f  m  : i  '  t'j  v  •  ;  r 

»  .>  j|Hf  $  i\i  ;  j.  W 
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I  Theorem  9  together  with  the  preceding  remarks  prove  the  final  theorem. 

Theorem  10.  (SA)  also  solves  (SP)  tinder  (MP).  □ 

It  is  easy  to  get  a  good  bound  on  the  running  time  of  the  combinatorial  part 
of  both  (PA)  and  (SA).  Let  v  be  the  number  of  non-zeros  in  A,  which  we  can 
l  assume  is  greater  than  ».  We  use  the  following  trick  to  reduce  the  running  time 

•  of  both  (PA)  and  (SA).  For  (SA)  as  well  as  (PA),  find  a  fixed  initial  maximum 

■  matching  on  A;  this  takes  0{mv )  operations.  Then,  when  finding  a  maximum 

?  matching  in  a  zero-section,  copy  over  the  part  of  the  fixed  matching  that  lies 

'f,  in  the  columns  of  the  zero-section  as  the  starting  matching;  this  copying  takes 

f  0{m?)  time.  Since  every  initially  unmatched  row  in  the  zero-section  matches 

to  some  column  outside  the  zero-section  in  the  fixed  matching,  the  number  of 
I  unmatched  rows  in  the  starting  solution  for  row  *’s  zero-section  can  be  at  most 

•j  the  number  of  non-zeros  in  row  t.  Thus  the  total  number  of  augmentations 

;■  needed  over  all  rows  is  0{v).  Since  each  augmentation  is  an  0{v)  operation,  we 

get  an  Ofv2)  overall  bound  for  the  combinatorics. 
j  The  time  needed  to  do  the  numerical  part  of  (PA)  and  (SA)  can  be  bounded 

■J  as  follows.  In  the  worst  case  we  will  have  to  solve  a  linear  system  like  (4)  of 

dimension  0{m)  for  each  one  of  m  rows.  Solving  one  such  system  is  bounded  by 
ij  0(ma),  so  the  numerical  part  is  bounded  by  0(m4)  overall.  However,  we  have 

assumed  that  A  is  sparse,  and  there  are  sparse  equations  routines  that  can  solve 
j  a  system  like  (4)  in  time  more  like  0(m3).  In  practice  we  expect  to  see  only  0(1) 

!  rows  whose  linear  systems  are  really  as  large  as  0(m);  most  linear  systems  will 

j  be  of  size  0(1).  Thus  under  favorable  circumstances  the  numerical  computations 

J  could  take  time  as  small  as  0(ma). 


4.  PreetieaNties 


Very  few  real-life  matrices  satisfy  (MP).  In  light  of  Theorem  5  we  cannot 
hope  to  actually  solve  (SP)  on  all  real  matrices,  but  we  can  try  to  apply  one  of 
our  algorithms  or  a  variant  as  an  'optimal*  heuristic.  Ideally,  when  we  apply 
our  "real"  algorithms  to  real,  full-rank  matrices,  they  would  be  guaranteed  to 
achieve  at  least  the  increase  in  sparsity  that  an  “ideal"  algorithm  would  achieve 
on  a  matrix  with  the  same  sparsity  pattern  that  did  satisfy  (MP). 

It  is  difficult  to  anticipate  unexpected  cancellation  with  real  matrices.  A 
parallel  type  of  algorithm  is  therefore  unsuitable,  as  it  has  to  proceed  without 
knowing  where  cancellation  takes  place.  On  the  other  hand,  a  sequential  type  of 
algorithm  can  take  stock  of  the  cancellation  that  arises  at  each  step.  However, 
guaranteeing  performance  becomes  more  subtle  in  the  presence  of  cancellation. 
Consider  the  full-rank  matrix 

(1  3  0  5  5\ 

2  1  4  0  0. 

0  3  0  5  5/ 
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Further  Question*  ud  Conclusion 


9 


I 


► 


Any  sequential  algorithm  ■will  pick  Ui  =  {  3  },  and  could  pick  G\  —  {  2  }.  This 
transformation  unexpectedly  seros  out  columns  4  and  5  of  row  1.  Thus  if  we 
naively  process  row  2  using  this  new  row  1,  we  will  choose  Ui  —  { 1 }.  But 
the  parallel  U\  —  0,  which  does  not  contain  U%  as  required  by  the  induction 
hypothesis  of  Theorem  8,  so  we  can  no  longer  guarantee  that  our  final  answer  will 
be  as  good  as  the  ideal.  (A  close  reading  of  that  proof  of  Theorem  8  will  reveal 
that  the  only  way  that  this  difficulty  can  arise  is  when  rank  A**  <  M{QRZ)  for 
some  sero-section;  in  the  second  zero- section  of  this  example,  rank  ARZ  =  1  < 
2  =  M{$rZ).) 

A  simple  trick  will  avoid  this  problem.  As  we  perform  (SA),  we  certainly 
know  at  each  step  where  we  expect  the  non-seros  to  occur  for  subsequent  steps. 
If  we  do  encounter  unexpected  cancellation,  we  can  merely  pretend  that  there 
is  still  a  non-zero  in  the  cancelled  position.  That  is,  subsequent  matchings  are 
performed  as  if  no  unexpected  cancellation  ever  took  place,  though  we  keep  track 
of  which  “non-zeros”  are  really  zeros.  Then  the  proof  of  Theorem  8  becomes 
valid  once  again,  and  the  modified  (SA)  is  now  guaranteed  to  produce  an  answer 
at  least  as  good  as  the  "ideal”  answer. 

We  now  make  some  remarks  about  implementing  (SA).  Linear  constraints 
are  usually  presented  as  a  mixture  of  equalities  and  inequalities.  If  these  con¬ 
straints  are  converted  to  the  form  (1)  by  adding  a  slack  variable  to  each  inequality 
row,  it  is  easy  to  see  that  there  is  always  a  maximum  cardinality  matching  in 
which  every  inequality  row  is  matched  to  its  slack  column.  It  is  also  easy  to 
see  that  such  rows  can  never  be  profitably  used  in  the  optimal  U  for  any  other 
row,  since  the  slack  variable  will  always  unavoidably  fill-in  its  column.  (In  fact, 
by  this  same  reasoning,  if  A  is  known  to  have  an  embedded  identity  matrix, 
then  A  must  already  be  optimally  sparse.)  Hence  (SA)  will  still  work  correctly 
if  we  merely  treat  inequality  rows  as  if  they  were  unmatched,  without  having 
to  explicitly  create  a  slack  variable  at  all.  This  phenomenon  implies  that  (SA) 
will  tend  to  find  better  solutions  for  systems  with  a  high  proportion  of  equality 
constraints. 


S.  Farther  Questions  end  Conclusion 

Trying  to  solve  the  Sparsity  Problem  as  described  in  this  paper  ruses  some 
interesting  questions.  From  an  applications  point  of  view,  the  chief  question  is: 
Does  (SA)  help  in  practice  or  not?  The  answer  to  tins  question  must  come  from 
empirical  experience  with  (SA)  on  various  problems.  We  have  implemented  a 
preliminary  version  of  (SA)  for  this  purpose;  our  results  so  far  are  encouraging, 
but  we  have  by  no  means  conclusively  demonstrated  the  usefulness  of  (SA).  We 
expect  to  report  our  computational  experience  with  (SA)  in  the  near  future. 

Although  it  is  necessary  to  keep  unexpected  zeros  as  phantom  non-zeros 
to  guarantee  the  performance  of  (SA)  on  real  matrices,  it  is  certainly  feasible 
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to  ran  the  algorithm  without  this  artifice.  Can  any  guarantees  be  made  in 
this  case?  Does  this  make  much  difference  in  practice?  Alternatively,  since 
"lucky*  cancellation  has  been  observed  in  nearly  all  real  examples  we  have  tried, 
is  there  some  efficient  heuristic  for  (SP)  that  can  take  advantage  of  this  and 
outperform  (SA),  perhaps  restricted  to  some  subset  of  interesting  problems  with 
special  structure?  Finally,  what  happens  when  we  try  to  apply  these  algorithms 
to  rank-deficient  matrices? 

We  shall  continue  our  research  on  (SP)  and  shall  try  to  answer  some  of  these 
questions  in  future  papers. 


Appendix. 

Proof,  (of  Theorem  1)  This  Theorem  and  its  proof  are  due  to  L.  Stockmeyer 
(1982).  See  Garey  and  Johnson  (1979)  for  the  definitions  of  the  concepts  used  in 
this  proof. 

The  problem  that  we  shall  reduce  to  (SP)  is 

Simple  Max  Cut:  Given  an  undirected  graph  Q  —  (V,  E),  partition  the  nodes  of 
Q  into  P  and  V  \  P  so  as  to  maximise 

({{«•,/}  €E  |  i€P,j€V\P}l. 

Let  n  =  |V|,  m  =  \E\,  let  A(§)  be  the  usual  (0,1)  node- arc  incidence  matrix 
of  Q,  and  let  At  be  the  n  X  2m  matrix  which  is  all  sero  except  for  row  i,  which 
is  half  +1  and  half  — 1.  Let  t  be  the  2m- vector  of  all  ones  and  let  /  be  the 
(2m(n  + 1)  +  l)-vector  of  ones.  Now  suppose  we  could  solve  (SP)  on  the  matrix 

A.  M  l  $} 

Clearly  some  row  of  the  optimal  T  for  B(Q),  which  we  may  assume  without  loss 
of  generality  is  the  first,  will  be  a  multiple  of  (1,  «i,  ca,...,  «»)  where  =  ±1  for 
all  t'  €  V.  (Because  of  the  sise  of  /  the  first  column  of  T  must  be  (1,0,0, . .  .,0) 
and  so  this  choice  for  the  first  row  of  T  causes  no  singularity  problems.)  Let 
P  =  {i\€{  =  +1}.  Then  the  number  of  non-seros  in  the  first  row  of  B($)  is 
clearly 

(2m(!»  +  l)  +  l)  +  na  +  (m-K{<1/}€£|<ei,,/€V\i>}|).  (7) 

But  since  (7)  is  minimised  by  the  optimal  T,  P  also  solves  the  Simple  Max  Cut 
Problem  for  Q.  □ 


RcfbnneM. 

Bunch,  J.  R.  and  D.  J.  Rose,  eds.,  “Sparse  Matrix  Computations,”  Academic 
Press  (New  York,  1976). 

Duff,  1.  S.,  ed.,  “Sparse  Matrices  and  their  Uses,”  Academic  Preu  (New  York, 
1980). 

Ford,  L.  R.  and  D.  R.  Fulkerson,  “Flows  in  Networks,”  Princeton  University 
Press  (Princeton,  1962). 

Garey,  M.  L.  and  D.  S.  Johnson,  “Computers  and  Intractibility,”  Freeman  (San 
Francisco,  1979). 

Lawler,  E.  L.,  “Combinatorial  Optimisation,”  Holt,  Rinehart  and  Winston  (New 
York  1976). 

Papadimitriou,  C.  H.  and  K.  Stieglits,  “Combinatorial  Optimisation:  Algorithms 
and  Complexity,"  Prentice  Hall  (Englewood  Cliffs,  1982). 

Stockmeyer,  L.  J.,  personal  communication  (1982). 


‘m-.V. 


REPORT  DOCUMENTATION  PAGE 


SOL  82-13 


4.  TITU  (m4  MWi) 


A  Fast  Algorithm  That  Makes  Matrices 
Optimally  Sparse 


ORIENT'S  CAT  A  LOO  NUMBER! 


».  TYPE  OP  REPORT  ft  PERIOO  COVKREO 


Technical  Report 


S.  HRPORWM  ORO.  REPORT  NUMBER 


AUTNO Mfe) 

Alan  3.  Hoffman 
S.  Thomas  McCormick 


M00014-75-C-0267. 

•m&it-'rt-c-oiio 


PERFONMINO  OROANIZATION  NMlK  AND  AD 

Department  of  Operations  Research  -  SOL 
Stanford  University 
Stanford,  CA  94305 


11.  CONTROLLINO  OFFICE  NAME  AND  ADDRESS 

Office  of  Naval  Research  -  Dept,  of  the  Navy 

800  N.  Quincy  Street 

Arlington,  VA  22217  _ 


Cmutnlllmg  OtBmm)  If.  SECURITY  CLASS,  (mi  M«  Mp«0 

UNCLASSIFIED 


NR-047- 143 


It.  REPORT  DATE 

September  1982 


IS.  DISTRIBUTION  STATEMENT  (ml  BU«  Rye  TO 


This  document  has  been  approved  for  public  release  and  sale; 
its  distribution  Is  unlimited. 


17.  DISTRIBUTION  STATEMENT  (mi  0»m 


IS.  KEY 


sparse  matrices 
bipartite  matching 
linear  constraints 


ABSTRACT 


SEE  REVERSE  SIDE 


7^3" 


SOL  82-13  A  FAST  ALGORITHM  THAT  MAKES  MATRICES  OPTIMALLY  SPARSE,  Alan  3 
Hoffman,  S.  Thoaas  McCormick,  (September  1982,  11  pp.) 


Under  a  non-degeneracy  assumption  on  the  non-zero  entries  of  a  given  sparse 
matrix,  a  polynomially-bounded  algorithm  is  presented  that  performs  row 
operations  on  the  given  matrix  which  reduce  it  to  a  sparsest  possible 
matrix  with  the  same  row  space.  For  each  row  of  the  matrix,  the  algorithm 
performs  a  maximum  cardinality  matching  on  the  bipartite  graph  associated 
with  a  submatrix  which  is  induced  by  that  row.  The  dual  of  the  optimal 
matching  then  specifies  the  row  operations  that  will  be  performed  on  that 
row.  We  also  describe  a  variant  algorithm  that  processes  the  matrix  in 
place,  thus  conserving  storage  and  time.  The  modifications  needed  to  apply 
the  algorithm  to  matrices  that  do  not  necessarily  satisfy  the 
non-degeneracy  assumption  are  also  described.  A  particularly  promising 
application  of  this  algorithm  is  in  the  reduction  of  linear  constraint 
matrices. 
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